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Abstract — In this paper, the Cramer-Rao lower bound (CRLB) 
for spatial correlation matrices is derived based on a rigorous 
model of the doubly selective fading channel for multiple-input 
multiple-output (MIMO) orthogonal frequency division multi- 
plexing (OFDM) systems. Adopting an orthogonal pilot pattern 
for multiple transmitting antennas and assuming independent 
samples along the time, the sample auto-correlation matrix of 
the channel response is complex Wishart distributed. Then, the 
maximum likelihood estimator (MLE) and the analytic expression 
of CRLB are derived by assuming that temporal and frequency 
correlations are known. Furthermore, lower bounds of total 
mean squared error (TMSE) and average mean squared error 
(AvgMSE) are deduced from CRLB for asymptotically infinite 
and finite signal-to-noise ratios (SNR's), respectively. According 
to the lower bound of AvgMSE, the amount of samples and 
the order of frequency selectivity show dominant impact on the 
accuracy of estimation. Besides, the number of pilot tones, SNR 
and normalized maximum Doppler spread together influence the 
effective order of frequency selectivity. Numerical simulations 
demonstrate the analytic results. 

Index Terms — CRLB, Spatial correlation matrices, Doubly 
selective fading channels, MIMO, OFDM, Complex Wishart. 

I. Introduction 

Due to the virtues of orthogonal frequency division mul- 
tiplexing (OFDM) for converting frequency selective fading 
channels into flat fading ones and of multiple-input multiple- 
output (MIMO) techniques for exploiting spatial diversity gain 
and/or enhancing the system capacity, many current systems 
combine them two together to achieve a better quality as well 
as a higher throughput [1]. 

For MIMO systems, the spatial correlation matrices play 
very important roles and are widely utilized, for example, to 
facilitate the transmitting precoding [2] [3], MMSE receiver 
[4] and multi-user strategies [5], However, since the true 
spatial correlation matrices are unknown in real applications, 
the sample correlation matrices have to be used in stead and 
are usually obtained through the channel estimation. 
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In this paper, we study the Cramer-Rao lower bound 
(CRLB) for the sample spatial correlation matrices of dou- 
bly selective fading channels for MIMO OFDM systems. 
Based on the a rigorous doubly fading channel model and 
assuming invariant pilot sequence along the time, the maxi- 
mum likelihood estimator (MLE) and the CRLB are derived. 
Then, the analytic expressions of lower bounds of the total 
mean squared error (TMSE) and average mean squared error 
(AvgMSE) are obtained for asymptotically infinite and finite 
signal-to-noise ratios (SNR's), respectively. Based on the lower 
bound of AvgMSE, several factors influencing the accuracy 
of estimation, including the amount of samples, the order 
of frequency selectivity, the normalized maximum Doppler 
spread, the number of pilot tones per antenna and SNR, are 
further analyzed. 

This paper is organized as follows. In Section [II] the MIMO 
OFDM system and channel model are introduced. Then, in 
Section [Till CRLB of the sample spatial correlation matrix is 
derived and further lower bounded to uncover the essential fac- 
tors. Numerical results appear in Section [IV] Finally, Section 
IVl concludes the paper. 

Notation: Lowercase and uppercase boldface letters denote 
column vectors and matrices, respectively. (•)*, (-) T , (-) H , 
(-)T, and || • \\p denote conjugate, transposition, conjugate 
transposition, Moore-Penrose pseudo-inverse and Frobenius 
norm, respectively. ® denotes the Kronecker product. E(-) 
represents expectation. [A]ij and [a]j denotes the (zj')-th 
element of A and the i-th element of a, respectively, diag(a) 
is a diagonal matrix by placing a on the diagonal. 

II. System Model 

Consider an MIMO OFDM system with a bandwidth of 
BW = 1/T Hz (T is the sampling period). N denotes the 
total number of tones, and a cyclic prefix (CP) of length L cp 
is inserted before each symbol to eliminate inter-block inter- 
ference. Thus the whole symbol duration is T s = (N + L cp )T. 
Then, tit transmitting antennas at the base station (BS) and 
tir receiving antennas at the mobile station (MS) are assumed, 
respectively. 



Between the i-th transmitting antenna and j-th receiving 
antenna, the complex baseband model of a linear time-variant 
mobile channel with it 3 '*' paths can be described by [6] 

h^(t,r)=Yl h!f A (t)5(r-r^T) (1) 
1=1 

where (t;)^ 4 ) G 1Z is the normalized non-sample-spaced 
delay of the l-th path, and ft; (t) is the corresponding 
complex amplitude with the power (of)"' 1 ). 

The following conditions are assumed to characterize the 
correlation property of the channel. 

1) Space Correlation: the stochastic MIMO radio channel 
model [7] is adopted; 

2) Frequency Correlation: the wide-sense stationary uncor- 
rected scattering (WSSUS) [6] is assumed; 

3) Time Correlation: the uniform scattering environment 
introduced by Clarke [8] is assumed. 

4) Scattering Function Separability: the channel has degen- 
eracy in all three dimensions [9]. 

Therefore, the complete spatial correlation matrix of the 
MIMO radio channel is given by [7] 

H s = H SiT <g> E StR (2) 

where S S , T & C" tX " t and E SvR G C nRXnR are the sym- 
metrical complex correlation matrices of antenna arrays of BS 
and MS, respectively. Then, the normalized time correlation 
function of any path is identical, i.e., [8] 

r t {At) = J Q (2irf d At) (3) 

where fd is the maximum Doppler spread, and Jq(-) is the 
zeroth order Bessel function of the first kind. In addition, 
LCM) = L, (n)^ = ti and (of)^ = of, hence, the 
frequency correlation matrix is 

R/=F r DFf (4) 

where D = diag(er, 2 ), l — 1, . . . , £ and F r G C NxL is the 
unbalanced Fourier transform matrix, defined as [F T ]^.; = 
e -j27rfcri/AT jy[ oreoveri jjjg power of the channel between each 
pair of transmitting and receiving antennas is normalized, i.e., 
tr(D) = 1. 

Assuming a sufficient CP, i.e., L cp > L, the discrete signal 
model in the frequency domain is written as 

tit 

yf (n) = Kf A ( n ) x f (") + n /' ) (") ( 5 ) 

where (n), (n), ^ (n) G Q Nxl are the n-th trans- 
mitted and received signal and additive white Gaussian noise 
(AWGN) vectors, respectively, and nf ' l \n) 6 C NxN is the 
channel transfer matrix between the i-th transmitting and j-th 
receiving antennas with the (k + v, fc)-th element as 

N-i L 

H^n)] = ± Yl E m)e-^ vm+kT '^ N 

J f5~\~u ,k 1\ _ 

m=0 1=1 

(6) 



where h\^ l \n,m) = h^' l \nT s + (L cp + m)T) is the sampled 
complex amplitude of the Z-th path, k and v denote frequency 
and Doppler spread, respectively. 

III. CRLB of Spatial Correlation Matrices 

Usually the correlation matrices of the channel response are 
obtained through the least squared (LS) channel estimation on 
pilot tones, that is, only pilot symbols are extracted and used 
to perform LS channel estimation. Therefore, 

h^ ) (n) = (X»(n))- 1 y? ) (n) (7) 

where X^(n) = diag(xi 1 ^ (n)) is a diagonal matrix consist- 
ing of pilot symbols. To alleviate the interference between 
multiple transmitting antennas, pilot symbols, i.e., x^ (n), 
i = 1, . . . , riT, are designed to be orthogonal, therefore 

(xf\n)) H xf >(n) = S(h - i 2 )\^f\n)\l (8) 

In this paper, we adopt the frequency division pilot pattern 
of which the pilot tones allocated to a certain transmitting 
antenna is exclusive of the others. To be more specifically, 
it is assumed that the pilot tone set of the i-th transmitting 
antenna is T p ^ = {i + kx9;k = 0, . . . , P — 1}, where P is the 
number of pilot tones of each antenna satisfying L<P, 9>riT 
and P9 < N. Thus, the i-th transmitting antenna transmits 

(i) 

non-zero pilots on T p while nulls on the rest. Besides, we 
assume that the non-zero elements of x^ (n) constitute a fixed 
vector x p G C Pxl , which is independent of n and i and of 
the normalized power so that X p X^ = Ip, where X p = 
diag(xp). Then, by discarding the zero elements, is further 
rewritten into 

h$>(n) = X; 1 H«.^)(n)x p + X^n^V) (9) 

where the subscript p denotes elements on pilot tones. 

According to the orthogonal pilot pattern, the instanta- 
neous channel impulse response (CIR) vector between the i- 
th transmitting antenna and the j-th receiving antenna cor- 
responding to the (i + kx9)-th sample of the n-th OFDM 
symbol can be denoted as hj % \n, i + kx9) = [h^ (n, i + 
kx6),...,h i i' i) (n,i + kx9)] T , k = 0, ...,P- 1. Then, the 
(j, i)-th CIR matrix for the n-th OFDM symbol is formed as 
H^ j) (n) = \b^ A (n,i),...,b^' i) (n,i + (P - 1) x <?)]. Ac- 
cording to the assumptions of WSSUS and uniform scattering, 
H[ j ' l \n) is complex normal, i.e., 

H^' l) (n) ~CAA L xp(0,O®D) (10) 

where fi G C PxP is a Toeplitz matrix, defined as 

[n} klM =J (27rf d (k 1 -k 2 )9T) (11) 

Then according to © and the pilot pattern, the (j, i)-th 
channel transfer matrix is H p ' l \n) — Ff'Hj^(n), where 
Fr is the submatrix of F r by drawing rows from Ip . Thus 

HW'fl(n) ~Ctf PxP (Q,n® [F«D(F«) H ]) (12) 



Assuming CIR is independent of the thermal noise, with (0 
and ( fT2b . we have 



MA) 
P-.ls 



(n) ~CM p $,tM) 



(13) 



where the covariance matrix T^''*J is defined as 

TO'.*) = (xf Ox^X^ F«D(F«) ff X p ) + CTl 2 j p (14) 

Apparently, T! 3 '' 1 ' is irrelevant of the indexes of the receiving 
antennas, which follows the assumption |4). According to the 
orthogonal pilot pattern, = F^**, where G C PxL 
is a submatrix of F T by drawing rows from If 1 = {k x 
6; k = 0, . . . , P — 1} and $ is a diagonal phase-twisted matrix 
with [3?]u = e~i 2lTTl / N . Then, the frequency auto-correlation 
matrix of is 

RW = F«D(FW) H = F^* i D(F^* l ) H = R p (15) 



where R p = F?'- ) D(FV U '')' H . Furthermore, to simplify the 
following derivation, we assume that CFR varies slightly 
within riT contiguous tones, so that the frequency cross- 



correlation matrix between the pilot sets Xp and 2 P 

R (u,*2) _ F (ii) D ( F (i2))H _ F (0)^ix-i 2D ^ F (0))J? ^ Rp 

(16) 

Then, the complete LS-estimated channel transfer matrix is 
constructed as 



r(ia) 



H Pi i s (n) 



v 



and it is complex normal, i.e., 



(17) 



' PtirXtit 



(18) 



where S G C PnRnTxPnRnT . According to the assumptions 
[D-0J, its ((ii — l)nR + ji, (*2 — l)nn + j2)-th submatrix is 
dl9l . shown at bottom of the next page. Then, with (15[ and 
(1 161 . X is expressed as 



S = H s (g (wA) + cr^I nr „ B p 



(20) 



where w = x^rix p and A = X^R p X p . 

When N t complete LS estimated channel transfer matrices 
are available, the sample auto-correlation matrix is formed as 



1 



Aft 

S = -^r vec(H p , is (n))vec(H p ^(7i)) ff 

* n— 1 



(21) 



where JVt is the number of samples. To derive the probability 
density function (PDF) of the sample auto-correlation matrix, 
we assume that samples are independent of each other, which 
may be a strict constraint. However, when the spacing between 
two contiguous pilot symbols is sufficiently large, the corre- 
lation between them is rather low, which alleviates the effect 
of model mismatch. Then, X has the complex central Wishart 



distribution with Nt degrees of freedom and covariance matrix 
S' = £/A^ [10], denoted as 



S~CW nTnsP (JV t) S') 



(22) 



and its PDF is 



J[ ' CT„ T „ H p(^ t )[det(S')]^ { ' 

where etr(-) = exp(tr(-)) and CT nTna p(N t ) is the complex 
multivariate gamma function, defined as 

n^n^P 

CT nTniiP {N t ) = n n T n R P(n T n R P-l)/2 JJ - k + 1) 

fc=l 

Then, according to i23\ . the likelihood function with respect 
to the parameter matrix H s can be written as 

L(S.) = trC-S'- 1 !:) + (N t - n T n R P) ln(det(E)) 

- \n{CT nTnR p(N t )) - N t ln(det(E')) (24) 

Therefore, the score function [11] is 

9L(B S ) 0vec(S') r dL(B s ) 
S co re (vec( as )) = = dyec( ^ ) 9yec(s0 (25) 

where the first multiplicative term on the right-hand side, 
according to ( 120) , is 



<9vec(£') i lu <9vec(E s ® A) J 



9vec(H s ) iVf 9vec(H s ) 
and the second term is 



(26) 



dL(S s ) 



<9vec(E') 



J - = vec[(S'- 1 ES'- 1 - Nil'- 1 ) 1 ] (27) 



Since 



vec(H s ® A) = K® [vec(H s ) <g> vec(A)] (28) 
where K a is 



K 



(29) 



where K (nR „ T)P G 7^(«««T)Px(n B nr)P is a transpose ma- 
trix satisfying vec(B T ) = K(„ H „ T )pvec(B), where B G 
C ni?nTxP . Hence, d26j is rewritten into 



<9vec(E 



AT 



avec(Hj = ^ [I (™^ Vec(A) ]K « (30) 

By letting the score function equal zero, we know that the 
maximum likelihood estimator (MLE) of E is E, so, with ( fT9l , 

its (fci,/c2)-th submatrix, denoted as {Ej-jj.^^ G C PxP , can 
be used to estimate [85]^.^ by 

[S s ] fcl , fc2 x (uA klM ) + a^h - k 2 )I P = m£ x p (31) 

where ki = (ii - l)np + j\, k 2 = (i% - l)np + j 2 , and 

A klM = XfR^)X p (32) 

Note that Ak x ,k 2 an d ^ assumed to be known for d3Tl ). 
that is, the frequency auto/cross-correlation matrices, the pilot 
sequence and the noise power are available. In order to solve 



[S s ]/si,fe 2 from (TJTI i, the singular matrices of fc a , denoted 
as XJk ly k 2 an d Vjfei ,fc 2 , is used to transform OTb into 



w[S.] fcl)fca A fcl)fc9 +^«(A!i-A a )Ii> = Ug ifa {E}jgfV fcl)fca 

where A^j.,, is a diagonal matrix with singular values 
of A^ £ 2 on the diagonal. Then, since rank^A^ _j. 2 ) = 
rank(Rp) = L, the MLE of [S s ]fe li fc 2 is 



MLE([E s ] fcl , fc2 ) = 

^ [Ug. fc2 {S}g x fc2 P) V fcl ^ 2 ] u - alSjk, - k 2 ) ^ 

2_^ Cl ' n i ( 33 ) 

1=1 



[A fel 



k 2 1,1 



where c;'s are normalized non-negative weight coefficients, 
i.e., q > and J2i=i c i = L 

Further, according to the score function, the Fisher Infor- 
mation matrix with respect to S s [11] is 



J(H S ) = E 



( <9L(E S ) \ / 9L(H 



\<9vec(H s ) J \(9vec(H 



H 



(34) 



with (|25]l(|30]t(|27]i, (fJUl is rewritten into ([35]), shown at the 
bottom of the next page, where B = E' _1 EE' _1 . Notice that 



B~CW N (N u -£' 



therefore 



E{vec[{B-N t T,'- 1 ) T ]vec[(B-N t 'E'- 1 ) T ] H } = Var[vec(B T )] 

(36) 

According to [12], (f36l > is 



Var[vec(B T )] = N t {H'' H ® £'~ T 
Then, with 07), J(H S ) is 

J(S«) = ^IVh-'t) 2 ®vec(A fcl , fc2 ) T ]K^(£ 

xK ®[I(n R n T ) 2 ® vec ( A fei,fc 2 )*] 

Therefore, the CRLB of H s is [11] 

CRLB(H S ) = J _1 (H S ) 



i-H 



(37) 

(38) 
(39) 



Now we consider the case that the SNR is asymptotically 
infinite, or, equivalently, the power of noise is zero. According 
to ( 1201 ). then, £' is reduced to 



S' = ^H s 



(40) 



Since A is rank deficient if L < P, E' 1 should be replaced 
by E' f . Then 



N? „ 



)A tff ®A tT (41) 



With dH]), (ED is rewritten into 

J(S.) = iY t [I (nj?nr)2 ®vec(A) T ][(H7 ff ®H; T ) 
®(At* ® At T )][I (nfl „ T)2 ® vec(A)*] 
= aAr t (H; H ®H s - T ) (42) 

where 

a = vec(A) T (A tH ® A tT )vec(A)* = L (43) 

From (l42l and (l43l , (1391 is rewritten into 

1 



CRLB(H S ) 



(44) 



Based on d44l i. a lower bound of TMSE of estimating H s is 
obtained, i.e., 

TMSE is (H s ) = tr[CRLB(H s )] = ( ™ )2 (45) 

And, accordingly, the lower bound of AvgMSE of H s is 

TMSE LB (H 5 ) 1 



AvgMSE iB (H s ) 



(n T n R ) 



(46) 



In real applications, the number of significant eigenvalues 
of A, denoted as L s , may be less than L. Since the amount of 
samples, N t , is finite, the insignificant eigenvalues are much 
less credible than significant ones. Therefore, only significant 
ones are used in d33l l. Besides, the weights of significant ones 
are considered to be equal. Hence, ( f33b is changed into 

MLE([H s ] fel , fe2 ) = 

^ [Uf^lSlifr^^M - fjWh - fca) (4?) 



£«[ A ki,ka] 



Correspondingly, the lower bound of average MSE of S is 
modified into 



AvgMSE LB (S a ) = 



1 



1 



(48) 



aN t L s N t 

In fact, L s represents the order of frequency selectivity, that 
is, the number of equivalent independent parallel transmission 
branches of the multipath channels for OFDM systems. 
When the SNR is finite, E' is approximated as 

(49) 



£' = ^E S ®(A + ^I P ) 

l\f to 



and, correspondingly, d39b is rewritten into 



CRLB(H S 



and d46i l is rewritten into 

AvgMSE iB (H s 



0N t 



(50) 



(51) 



&}V-lL+h,(i 2 -l)n R+j2 = [*s}(n-l)nn+n,(iz-l)nn +h x (u>X?R^X p ) + a^h - i2)5(h j 2 )l P (19) 



where 

p = vec(A) T [(A 

L 1 



^l P ) 



-H 



(A+^Ip) 

It? 



~ T ]vec(A)* 
(52) 



where p\ 
since 



and A; is the Z-th eigenvalue of A. Moreover, 



\l x R(H) (53) 



where ||x p ||2 is the power of pilot symbol, and R(J~2) is the 
Rayleigh quotient of J~i [13]. Due to the normalized power, 
Il x pll2 = P- Besides, it is straightforward that R(fi) < 
Amoi(^), where X ma x(^) denotes the maximum eigenvalue 
of 1~2. According to [14], when 6fdT s < 0.35, X max (Cl) can 
be well approximated by 



« PJ (2<Kc6f d T s ) 
where c = 0.35. Therefore, j3 is upper bounded by 

1 



L 



1=1 



(P?J (2ircef d T s ) Pl ) 
Then, ( BTT i is further lower bounded by 
AvgMSE LB (H s ) = 



A, 



PmaxNt 



(54) 



(55) 



(56) 



Note that ujpi can be regarded as the effective SNR on the Z-th 
subchannel. When uipi is too small, say, below dB, it should 
not be used in MLE (|47| |. which would reduce the effective 
order of frequency selectivity. According to (l55l l. therefore, 
the number of pilot tones, SNR and maximum Doppler spread 
together influence the effective order of frequency selectivity 
and, further, the accuracy of estimation. 

IV. Numerical Results 
The OFDM system in simulations is of BW = 1.25 MHz 
(T = l/BW = 800 ns), N = 128, and L cp = 16. Two 3GPP 
E-UTRA channel models are adopted: Extended Vehicular 
A model (EVA) and Extended Typical Urban model (ETU) 
[15]. The excess tap delay of EVA is [0, 30, 150, 310, 370, 
710, 1090, 1730, 2510] ns, and its relative power is [0.0, 
-1.5, -1.4, -3.6, -0.6, -9.1, -7.0, -12.0, -16.9] dB. 
For ETU, they are [0, 50, 120, 200, 230, 500, 1600, 2300, 
5000] ns and [-1.0, -1.0, -1.0, 0.0, 0.0, 0.0, -3.0, -5.0, 
—7.0] dB, respectively. The classic Doppler spectrum, i.e., 
Jakes' spectrum [6], is applied to generate the Rayleigh fading 
channel. The MIMO configuration is 4 x 4, and the correlation 
matrices of transmitting and receiving antennas are shown at 
the bottom of the next page, respectively [7]. Besides, the 
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Fig. 1 . Comparison of analytic results of d48t and numerical results for EVA 
and ETU channels when SNR — * oo and = 100Hz. 



number of pilot tones per transmitting antenna is 16, the pilot 
spacing is 8, and the pilot symbols are separated far enough 
from others to reduce the correlation among them. 

In FiglU we compare the analytic results d48b and the 
numerical results over a range of iV t 's for EVA and ETU 
channels, respectively, when the SNR is asymptotically infinite 
and fd = 100Hz. The pilot sequences are QPSK modulated 
and randomly chosen. Besides, different L s 's are tested to 
demonstrate that d48b is a tight lower bound of d47l ). Appar- 
ently, the analytic results meet the numerical ones quite well. 

In FigJU different SNR's and maximum Doppler spreads 
are tested to demonstrate that (TS6b is tight lower bound of (1471 ) 
when the pilot sequence, x p , is properly selected, specifically, 
the eigenvector of O associated with the maximum eigenvalue. 
Since the number of pilot tones, SNR and normalized max- 
imum Doppler spread together influence the effective order 
of frequency selectivity, the value of L s varies for different 
cases, which is reflected by varying AvgMSE. Further, it is 
obvious that SNR has a more significant impact on L s than 
the maximum Doppler spread, when the number of pilot tones 
are large enough, e.g., over 16. 



J(S S ) = -J [l {nRnT y ® vec(A) T ]K^{vec[(B - AT t S'- 1 ) T ]vec[(B - V t S'- 1 ) T ] H }K 8 [I (llHnT)2 ® vec(A)*] (35) 
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Fig. 2. Comparison of analytic results of )56t and numerical results for EVA and ETU channels under different SNR's and Doppler's. 




V. Conclusion 

In this paper, we derive the CRLB of spatial correlation ma- 
trices based on a rigorous model of the doubly selective fading 
channel for MIMO OFDM systems. With several necessary 
assumptions, the sample auto-correlation matrix of the channel 
response is complex Wishart distributed. Then, the maximum 
likelihood estimator is obtained, and the analytic expressions 
of CRLB as well as lower bounds of TMSE and AvgMSE 
are deduced for asymptotically infinite and finite SNR's, 
respectively. According to the lower bound of AvgMSE, the 
amount of samples and the order of frequency selectivity 
influence the accuracy of estimation dominantly. Besides, the 
number of pilot tones, SNR and maximum Doppler spread 
have effects on the effective order of frequency selectivity. 
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